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NOMENCLATURE 


h 

K 

P 

Q 

q 

T 


T 

t 


Specific heat 

Heat transfer coefficient 
Thermal conductivity 
Power 

Energy input 
Rate of heat flux 
Temperature 

Fluid temperature for convective 
boundary conditions . 

Ambient temperature 

Time 


j/kg K 

w/m 2 K 

w/mK 

w 

w 

w/m 2 

K 

K 

K 

sec 


Greek Symbols . 
a Thermal diffusivity 

P Mass density 


V 
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Vector gradient operator 

A constant which denotes the type 
of finite difference scheme used. 


M 2 /hr 

kg/m 3 


M 


-1 



Subscripts . 


av average 

e Element index 

i Node index 

0 Node index 

m Melting 

s Solid 

Superscripts . 

. Derivative with respect to time . 



ABSTRACT 


A thermal model has been developed for 
predicting the material removal rate during Electro 
Discharge Machining (EDM) . The transient heat conduction 
equation for the work piece which accounts for the heat 
absorption due to melting, has been solved by the 
finite element method. Gaussian heat distribution of 
power within a spark has been considered and the 
periodical updating of the shape of the work piece surface 
has been performed for the first time, unlike in earlier 
theoretical works. Computations have been performed for a 
single spark in the form of a pulse and a multi spark 
situation with a continuous power input. The effect of 
different variables on the material removal have been 
calculated and plotted. The observed results seem to match 
fairly well with the available experimental results. 



CHAPTER -I 


INTRODUCTION AND LITERATURE SURVEY 


1 .1 GENERAL BACKGROUND : 

There has been a rapid growth in the development 
of harder and difficult -to -machine metals and alloys 
during the last two decades . Conventional machining is 
uneconomical for such materials and the degree of 
accuracy and surface finish attainable is poor. Merchant (I960) 
emphasized the need for the development of newer concepts 
in the machining of metals in order to meet the challenges 
posed by modern engineering materials. By utilizing the 
results of basic and applied research, newer machining 
processes, often called as unconventional machining methods 
have been developed . 

Electric Discharge Machining (EDM) is the process 
of removal of material from the electrically conducting 
materials by electrical discharges between two electrodes 
i.e., anodic work piece and cathodic tool. The inter- 
electrode gap is filled with a dielectric medium. The aim 
of the process is to achieve controlled removal of material 
from the work piece . 
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Fig. 1 ,1 Classification of EDM processes [7 ] 










EDM has found wide applications in the machining 
of hard metals and alloys which cannot be machined easily 
by conventional methods. It thus plays a major role in the 
manufacture of dies, tools etc., made of tungsten carbides, 
stellites or hard steels . Alloys used in the aeronautical 
industry, for example, hast alloy, nimonic etc. can be 
machined conveniently by this process . This process has an 
added advantage of being capable of machining complicated 
components from hard materials to a high degree of accuracy 
and surface finish . 

It is important to predict the obtainable work 
profile from the given machining conditions and specified 
tool shape . This can be done by knowing the exact mechanism 
of metal removal and a computational procedure for the same . 

The metal removal in EDM is considered to be due to 
melting and vaporization of anode material . Many attempts 
[1,2, 3, 4, 5, 6 ] have been made by researchers, to theoretically 
compute the MRR. But, discrepancy between the analytical and 
experimental results has been observed as discussed below. 

1 .2 REVIEW OF PREVIOUS WORK : 

Many researchers [1 ,2, 3, 4,5, 6] have contributed to 
the study of the EDM process. Most of these researchers 
considered a single spark in their study and then generalized 
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"the analysis for the whole process. Before proceeding 
to discuss our own work, we shall briefly present the 
work done by earlier researchers in this area . 

C.C. Marty [1 ] demonstrated that for a given heat 
input and a maximum surface temperature , the crater depth 
and radius of the heat source are correlated by an 
universal function. His experimental data show that for many 
materials, the surface temperature is of the same magnitude 
as the boiling point. Also, the surface temperature in the 
centre of the spot is practically independent of heat source 
diameter. 

Y. Mukoyama [ 2 ] applied the theory of thermal 
conduction and calculated the relation between the discharge 
energy, discharge duration, material properties of the 
electrode and the melting depth at a discharge point. Material 
removal is considered to be by thermal phenomena and 
temperature distribution is based on heat conduction 
equation in cylindrical co-ordinates. The differential 
equation for temperature distribution is 


3T 

3T 


ai 


■ 3 2 T 

— - 

3r^ 


+ 1 
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3T 

3r- 


+ i!$) 


dZ 
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The difference between heat supplied and heat taken 
up due to latent heat of melting is considered to be the 
actual heat input at the discharge point. It has been observed 
that crater depth increases with increase in discharge duration 
upto a certain value and then drops. This Is true for all 
electrode materials . He showed that melting depth is 
proportional to material removal . Tool wear quantity is 
attributed to wear . Material removal efficiency is about 
5-15 '/. (15-60 ’/.for graphite tool). The analysis is simpli- 
fied due to the assumptions like, heat source size is constant 
regardless of discharge condition and the physical constants 
of electrodes taken are those at normal conditions. 

i 

Jilani and Pandey [3] examined the effect of plasma 
channel growth on metal removal and found a good correlation 
between theoretical and experimental results . They proved 
that in case of steel for a fixed pulse energy, volumetric 
erosion is not dependent on the current pulse form, but the 
wear ratio is influenced by the current waveform. They 
concluded that the theoretical depth to diameter ratio of a 
crater lies in the range of 3.5 - 34 ’/.and experimental 
3-30 */. . In order to arrive at these conclusions, they 
considered a disc shaped heat source and material removal by 
thermal erosion due to sparking. The governing equation for 
temperature is 



R 2 a 


1 2 


T(r,z,t) - 

krc ~~ 
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d t exp( -z 2 / 4a t ) exp ■_ -r^(4at+R 2 ) , 
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For plasma channel growth, it has been assumed that 
for a spark, the surface temperature is constant. This 
model gives comparatively good results for MRR, since 
plasma channel growth with time and adapted thermal diffusivity 
(a) - to account for material removal at boiling point - are 
taken into consideration. Due to arbitrarily chosen spark 
diameter, the predicted MRR values are higher than the 
experimental ones . 

Dharmadnikari and Sharma [ 4 1 presented a thermal 
erosion model in which a finite number of heat sources 
are considered to determine the material removal rate for a 
single pulse. It has been considered that the heat source 
which is initially point shaped to be spreading progressively 
over the crater surface . The energy is assumed to be 
distributed equally among the sparks and craters formed are 
hemispherical in shape . Material removal per pulse is given 
as equal to 

(8n a 3/Z t p ) (K 3 /Kp) P av /J 
MRR = (8« a l /Z p) (42/Kp)P av /J 
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where , 


P = 

V 


t (pulse time) 
duty factor = —=■ 


t (cycle time) 

V 


3 P^/ZJttp) 0 -’ 


= a constant for a given value of 


A Erden and B Kaftanoglu [5 ] presented two new 
models. one is based on the uniformly distributed point 
sources assumption. The other is an approximate, but general 
model, involving time variation of the heat source and its 
radius. That is temperature is considered to be a function 
of heat source density and source radius at a given time. The 
heat source density and radius are considered to be discrete 
step functions and based on this continuously variable heat 
source is simplified to an energy source of constant 
magnitude . 

The governing equation (for constant power) is 
taken to be 

T(r ’ t; - erf c ] 

Out of the many factors affecting the material removal, 
the authors [5] considered only pulse time power and electrode 
material. Agreement with the experimental data is shown to 
be satisfactory for all cases . 
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R . Snoeys and F . Van Dyck [6 ] considered melting 
and vaporization of metal to be dominating the material 
removal process in EDM. They considered that the portions of 
the metal where the temperature is larger than or equal to 
the melting value . Constitute the material which is removed 
by the EDM process. Their heat transfer model is treated in 
cylindrical coordinates with no internal heat generation and 
no mass flow. In order to take into account the latent heat, 
thermal diffusivity is taken as a function of melting heat, 
time and temperature the governing equation for temperature 
distribution is 

T ( r,z,t) - £ J Q (X n r) (>*[ erf(X n ,/It + ^_>-1 ] 

+ e n ' [erf (X at - ~=)+1 ] 1 
n 2 Vat J 

1 

where a = K/C - adapted thermal diffusivity 

Because no vaporization is accounted for^350 to 500 */. 
difference is observed between calculated and theoretical 
upper limit (at 100 ’/.efficiency). 

From above literature survey, it is evident that 
these models do not account for multi sparks, moving work 
surfaces ( in actual machining, tool and work surfaces are 
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moving with time), convection of heat to the dielectric 
and environment, effect of contamination on the dielectric 
properties etc . 

1 .3 OBJECTIVES OF PRESENT WORK : 

The main objective of the present work is to 
calculate the material removal rate during electro-Discharge 
machining . 

‘Till now, researchers have considered the same 
surface for successive sparks thereby setting a condition 
that the problem is being repeated for each spark, literally 
with constant material removal per spark. In the present 
analysis, the surface is being changed (updated) between 
every two successive sparks by considering the melting 
isotherm as the new surface. Due to this, the material 
removal per spark will vary . 

To achieve this objective, the following steps are 
being proposed. 

1 . To develop a suitable heat transfer model since the 
material removal is a thermal process . 

2. To apply the finite element analysis for the above 

problem since numerical methods give fairly good 
results and also because the geometry of the problem 
becomes very complicated after partial machining. 
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3. T 0 find out the amount of material removal for 

a specified time, applying continuous sparking of 
known power with Gaussian distribution. 
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CHAPTER- II 


MATHEMATICAL ANALYSIS OF THE PROBLEM 


2.1 MECHANISM OF ELM : 

Any form of electric discharge is accompanied by 
some wear of the electrode material. The Electro Discharge 
Machining (EDM) is based on this spark erosion effect. 

During EDM, the electric discharge occurs through a gas, 
liquid or solid, called the dielectric. The necessary 
condition for a discharge is the ionization of the dielectric 
by which the molecules of the dielectric are split up into 
positively charged ions and negatively charged electrons. 

The following processes take place when a discharge 
occurs between the EDM tool (cathode) and the work piece 
(anode), through a dielectric medium. When a sufficiently 
high voltage is applied across the electrodes, the intensity 
of the electric field between them builds up. At some 
predetermined value of the voltage, the individual electrons 
break loose from the surface of the cathode and are impelled 
towards the anode under the influence of electrical field 
forces. While moving in the inter electrode space, the 
electrons collide with the neutral molecules of the dielectric, 
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detaching electrons from them and causing ionization. 

After some time, the ionization becomes such that a narrow 
channel of continuous conductivity is formed. Then there 
is a considerable flow of electrons along the channel to 
the anode, resulting in a momentary current impulse or 
discharge. The liberation of a large pulse of energy leads 
to the generation of high temperature (between 8000° and 
12,000°C), causing fusion or partial vaporization of the 
metal and the dielectric fluid at the poinu of discharge. 

The metal, in the form of liquid drops, is dispersed into the 
space surrounding the electrodes by the explosive pressure 
of the gaseous products in the discharge. This results in the 
formation of a tiny crater at the point of discharge in the 
work piece ' 1 ] . 

Comparatively less metal is eroded from the tool 
(cathode) than the work piece (anode) due to two reasons. 

1 . The momentum with which positive ions strike 
the cathode surface is much less than the 
momentum with which the electrons -stream impinges 
on the anode surface . 

2. A compressive force is generated on the cathode 
surface by the spark which helps reduce the tool 


wear . 
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A simple EDM set up is shown in Fig. 2.1. 

Typically the circuit consists of resistor and capacitor 
elements which draw electrical energy from a DC source. 

Due to the continuous. oscillatory response of the R .C . 
circuit, the capacitor is charged upto a threshold voltage, 
at which it discharges electrical energy in the form of a 
pulse through the circuit. The frequency and shape of the 
electrical pulse can be adjusted as desired by choosing 
proper values of resistance and capacitance for the circuit 
elements. 'Thus, a short pulse ( of the order of micro seconds) 
of very high energy intensity is generated which results 
in a spark discharge across the dielectric channel. 

2 .2 CONSIDERATIONS INVOLVED IN THE THEORETICAL 
ANALYSIS OF EDM ; 

The method for calculating the material removal 
rate due to the impact of an electric pulse, is based upon 
the assumption that the amount of molten metal right at the 
end of the an impulse is removed completely. The shock wave 
generated In the dielectric by the spark, and the rapidly 
expanding metallic vapour bubble at work surface, aid in 
clearing up the molten metal. In addition to these phenomenon, 
the dielectric fluid is often forced to flow through the 
inter -electrode gap (ieg) for effective material removal. 

Thus, the amount of molten metal can be determined by the 
volume limited by the iso-temperature plane of the melting 
temperature at the end of the pulse and the exposed surface of 
the electrode metal at the beginning of the spark. 
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During a typical electro discharge machining 
operation, for each electrical pulse, many sparks originate 
simultaneously from the tool electrode and they strike at 
the work-surface. Dharmad’nikari and Sarma " 4 ■ have found 
that the number of simultaneously produced sparks at any time 
is of the order of 15. However, it is the function of 
surface roughness of tool and work. The diameter of the 
spark itself is of the order 50 microns, while the heat 
affected area for each individual spark is about a hundred 
times larger than the spark area [ 4 . . The locations am which 

the sparks strike on the work-surface and the energy carried 
by each spark vary randomly and a precise estimation of these 
quantities is not possible. Also, after partial machining, 
the surface of the work-piece is filled with many complex- 
shaped craters. In view of these difficulties the considera- 
tion of all the sparks and their interaction with the work- 
metal is a very difficult proposition for analytical or 
numerical modelling of the EDM process. Many authors [ 1 , 2 , 3, 4, 5,6 
have considered a single, cylindrical heat source of uniform 
intensity or many point sources of equal power which are 
uniformly distributed over the area under the tool electrode. 
Mathematical analyses have also been performed with assumption 
that the location of heat source (or heat sources) is at 
the level of the original unmachined surface . 
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In the present work, a finite element analysis 
of the EDM process is performed by considering a given 
number (N) of sparks to strike at the work-surface. Due 
to the high pulse frequency, each spark is taken to be a 
continuous heat source. It is also assumed that the sparks 
strike at the same locations for all time and each spark 
is a cylindrical heat source of equal power. The intensity 
variation within the spark has been assumed to follow a 
Gaussian distribution. Considering a typical spark which is 
surrounded by other sparks , it can be taken that only a 
limited region is affected by the spark. The analysis, therefore, 
can be simplified by assuming the boundary of the heat 
affected zone of a particular spark to be loosing no heat 
to the surrounding zones, while calculating the rate of 
machining due to a single spark. Approximately, the total 
machining rate can be calculated by multiplying the individual 
machining rate by the number of sparks, N. 

2.3 HEAT TRANSFER MODEL FOR THE HEAT AFFECTED ZONE 
OF A SINGLE SPARK : 

In the present study each spark is considered to be 
cylindrical in shape as represented in Fig. 2.2. 

Considering the heating due to every spark to be 
axisymmetric in its neighbourhood, the simplification can be 
made that there is no temperature variation with respect to 


the direction 9 . 
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work piece 


Fig. 2.2 Cylindrical heat source 
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3T 

3 © 


= 0 


( 2 . 1 ) 


Under the above conditions, the governing heat 
conduction equation within the work material is given in 
cylindrical coordinates as 


P c T 


3T 

at 


= k 


a 

3r 




3 2 T 

7? 


( 2 . 2 ) 


where, p is the density, C is the specific heat and k is 

Jr 

the thermal conductivity of the work material . 

In order to account for the latent heat of melting 
and evaporation, the specific heat of the electrode 
material is taken to be a function of temperature as shown 
in Fig. 2.3. With this the governing equation (2.2) can be 
written as 


3T 


3T> 


0 V T > If - k [ fr Sf> + 


3 2 T 

3z 2 


] (2.3) 


Average values of different properties are used 
as their variations are not large in the range of temperatures 
of interest in this problem. 

The following boundary and initial conditions are 
to be satisfied, for the heat affected zone of a single 
spark, as shown in Fig. 2.4. 
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Temperature 


Cp g = Specific heat of solid 
Cpi = Specific heat of liquid 

T m = Melting temperature 
Fig. 2.3 Specific heat variation with temperature 




S, , Insulated boundary 

Heat flux boundary 
Convective boundary 


Fig. 2.4 Boundary conditions on the heat 
affected zone (HAZ). 
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The initial condition for temperature for the 
entire solution domain at the start of the machining 
operation is assumed to be the room temperature, T q . Thus 

T(r,z,o) = T q (2.4a) 

The boundary conditions are : 

(i) &11 time, the bottom boundary 3^ of the heat 

affected zone is considered to be insulated in view 
of the spark being surrounded by neighbouring 
sparks. Therefore, 


3T 

3n 


0 on 


(2.4b) 


where n is the coordinate normal to . 

(ii) Since the heating process is taken as axi -symmetric 
about the spark axis, 

~ = 0 at r = 0 (2.4c) 


(iii) On the top surface, heat flux is considered only 
upto a radius equal to that of the spark radius, 
Resurface S^) . For the rest of the portion (Surface S^), 
convective heat loss to the dielectric is assumed. 


- k = h (T-T f ) if r > R 
= Q(r) if r ^ R 


(2.4d) 
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where h is "the heat transfer coefficient between the 

work surface and the dielectric and Q(r) is the heat flux due 

to the spark. 

2.4 CALCULATION OF THE HEAT FLUX DUB TO SPARK : 

It is necessary to prescribe the distribution of 
energy intensity across the spark for the heat transfer 
model described above, while solving for the temperature in 
the heat affected zone . Many researchers have considered the 
spark energy to be uniformly distributed across the spark 
for mathematical modelling. However, this may not correspond 
to the actual case since the spark intensity depends on 
the surface irregularities (surface roughness). Even for a 
visibly smooth surface, some surface roughness will prevail 
which has a predominent effect since the spark itself has 
dimensions of the order of microns . Apart from determining 
the variation of energy intensity across the spark, it is 
also necessary to relate the overall parameters of the EDM 
process to the total power of a spark. This power per spark 
can be calculated as 

P s - -^SSSL (2 .5) 

where V is the average voltage, I is the average current, 
cos <$ is the duty factor and N is the number of sparks . 

It is also possible to derive this quantity by integrating 
the instantaneous discharge power over the pulse time . 
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Thus, 

P S - T- = B h * d at (2-6) 

d o 

where is the total energy supplied per pulse, t^ is 
the discharge time and q(t) is the instantaneous power. 

In the present work, a Gaussian distribution (Fig. 2. 5) 
of the energy intensity is considered within the spark. 

If q Q is the maximum intensity at the axis of the spark and 
R is the spark radius, for the Gaussian distribution, we get, 

q(r ) = q Q exp {- 2(|- ) 2 } (2.7) 

In terms of the total power per spark, the above 
expression may be rewritten as 

q(r) = exp { - 2 (|) 2 } (2.8) 

2.5 SOLUTION APPROACH TO THE SPARK EROSION PROCESS : 

Earlier authors [1,2, 3, A, 5, 6] have applied the 
mathematical solutions for the heat conduction problem in 
cylindrical geometry for both single and multiple spark 
models. These solutions have a major disadvantage that they 
consider the heat sources to be placed at the level of the 
unmachined surface for all times and do not consider the 
change in geometry after partial machining. The final 
analytical solution for a single heat source situation is 

provided in the form of an infinite series given by (See Ref .8) 
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R 



Fig. 2.5 Gaussian distribution of heat source 
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T(r,z,t) - T q = 


U 

£ n 


X n z 


n=1 ^r J o (X n r) L e :erf(X n ./at 


-X 2 

+ ^ -1 ( + e n { erf (X -/at — )+ 1 

2.'iaX 1 n 


( 2 . 9 ) 

A similar analytical solution can be obtained for 
a multispark situation also. The rate of machining in these 
models is predicted by plotting the melting isotherm as a 
function of time. 

In this work, the finite element approach has been 
applied to obtain more accurate predictions for MRR than the 
ones available at present. The major contribution of the work 
is the consideration of the change in shape of the work- 
surface after partial machining . After each computational 
machining cycle (which coincides with the time step used), 
the work surface is updated to the predicted shape of the 
melting isotherm of that time step. The entire finite 
element mesh is formulated and these calculations are carried 
out upto the desired total machining time. The details of 
the FEM calculations are presented in the next chapter. 
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2.6 LIMITATIONS OF THE PRESENT MODEL : 

• In actual practice, the EDM cycle consists of 

OFF time and ON (pulse) time. But in the present 
model the cycle is having only ON time thereby 
making it a continuous material removal process . 

2. The EDM, working cycle times are of the order of 
about 100b seconds. If it is to be taken into account, 
enormous computational time would be required to 
obtain machining for about a minute or so . In order 

to reduce the computer time, larger cycle times 
are taken. 

3 . Heat supplied to the evaporation of work material 
and dielectric are not considered. 

4. A prescribed number of sparks (twenty in all 
calculations), are assumed to strike between the 
interelectrode gap. In reality, this should be a 
function of surface roughness and various other 


factors 
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CHAPTER -III 

SOLUTION PROCEDURE 3Y FINITE ELEMENT METHOD 
3.1 INTRODUCTION : 

The Finite Element Method (FEM) is one of the 
numerical techniques for obtaining approximate solutions to 
a wide range of engineering problems . The FE method is widely 
used to solve a set of partial differential equations because 
it possesses certain advantages over the other numerical 
methods. For instance, if we compare the method with finite 
difference technique, the latter has the disadvantages such 
as - inaccuracy in approximating, the derivatives of the 
unknown variable and difficulty in handling complex and irregular 
geometries. General purpose software can be developed for 
FEM problems for a varied range of engineering and mathematical 
applications with different boundary conditions . Due to the 
above mentioned reasons, the present study is being treated 
by FEM technique for solving the heat transfer equations. 

In the Finite Element Method, the solution domain is 
conveniently divided into a finite, small interconnected 
sub -regions called elements. Solutions are formulated for 
each element individually and then assembled to form the 
solution of the entire problem. Since, we have choice over 
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the shape and grouping of elements, highly complex geometri- 
cal shapes can be dealt with. In FEM, the boundary contri- 
butions do not appear into the equations for the individual 
elements and the field variables need not be changed when 
the boundary conditions are changed. Also, a good approxi- 
mation to the boundary shape can be obtained due to the 
great flexibility In choosing the shape of the elements . 
Complicated material property variations can be incorporated 
into the system with ease . The systematic procedure for the 
derivation of the weighting functions makes it a powerful 
tool for dealing with non-linear and time -dependent field 
variables . 

3.2 DESCRIPTION OF THE FEM PROCEDURE ; 

The first step in FEM formulation is the conversion 
of the 3et of governing equations into an equivalent system 
of integral equations . Then the solution domain is divided 
into several small elements of chosen shape and size. The 
unknown field variables are expressed interms of assumed 
interpolation functions and the values of the field variables 
at specified points called nodes with each element. For the 
finite element representation of a problem, the nodal values 
of the field variables become the new unknowns . The integral 
governing equations of the problem are evaluated within each 
element using the interpolation functions for the element. 
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All such elemental integrals are assembled to obtain the final 
matrix equations . The unknown nodal variables can be 
calculated by solving the assembled matrix equations for the 
entire solution domain. 

Of the procedures available for converting the 
governing partial differential equations into equivalent 
integral equations ,the most popular ones are the variational 
procedure and the Galerkin'* s method. The variational approach 
relies on the extremization of a functional or a set of 
functionals, subject to the given boundary conditions. However, 
functionals may not exist for the problem at hand and the 
variational approach can not be applied to such problems . 

The Galerkin* s method or method of weighted residuals is more 
versatile in the sense that it can be extended even to highly 
non-linear problems where ho functional is available. In the 
present study the Galerkin’s approach is used because of 
its general applicability. 

3.3 GALERKIN METHOD : 

The Galerkin technique is one among the class of 
techniques collectively known as the Method of Weighted 
Residuals. In this approach, the governing differential 
equations are converted to integral form by minimizing the 
residue which arises due to the approximate solution. 



Iv/o basic seeps are involved in the application 
of the method of weighted residuals » the first stet is to 
assume the general functional behaviour of the dependent 
field variable in some way so as to approximately satisfy 
the given differential equation and boundary conditions . 
Substitution of this approximation into the original differen- 
tial equation and boundary conditions then results in some 
error at all locations which is known as the residual. This 
residual is required to vanish in some average sense over 
the entire solution domain. 

The second step in the method of weighted residuals 
is to solve the equations resulting from the first step and 
thereby specialize the general functional form to a 
particular function, which then becomes the approximate 
solution sought. 

Consider, 

/(GO - 0 (3.1) 

where 0 is the exact solution to the differential equation. 

An approximate solution 0 will leave behind a residue, R, 
such that 

/(/) = R (3.2) 

In Galerkin’s method, the residue weighted appro- 
priately by suitable weighting function, is minimized over 



"the whole solution domain u in an integral sense . In 
mathematical form, this will yield : 

; W R dv = s W H{(f) dv = 0 (3.3) 

D D x 

for each i = 1,2, m 

where VL are the weighting functions. The approximation 
0 may be written as a piecewise continuous profile for 
each element through the expression. 

= I N ^ (3.4) 

i=1 1 1 

where m Is the number of nodes per element, N^ e ^are the 
interpolation functions or shape functions and are 

the nodal values of within the element (e). In Galerkin’s 

procedure, the weighting functions are chosen to be the 

(e) 

same as the interpolation function N> J of each element. 

It is often advantageous to apply integration by 
parts to the residue minimization equation (3.3). By this 
technique, expressions containing lower-order derivatives 
are obtained which facilitate the use of approximating 
functions with lower-order inter -element continuity. 

Another advantage of integration by parts is that this 
offers a convenient way to introduce the natural boundary 
conditions of the problem. 
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3.4 APPLICATION OF PEM PROCEDURE TO THERMAL EROSION IN EDM ; 
3.4.1 Solution of two-dimensional heat conduction equation: 




where n is the total number of unknowns in the problem. 

The temperature in equation (3.7) can be expanded in terms 
of its nodal values and the interpolation functions in 
the form 

T = £ N. T. (3.8) 

3*1 J J 



2 2 


where are the interpolation functions. In Galerkln’s 
method, the weighting functions wh of equation (3.7) are 
taken to be the same as the inf erpolation functions * 
Therefore, 

1 1 N i l p C p H - K(^| + ^|) | dx dy =0 

3 x dy 

(3.9) 

Substituting equation (3.8) in equation (3,9) and applying 
integration by parts, 

3N. 3N. 3N. 3N- 

£ K//( 3ir 3t I+ 3T sf) dxd y] [Tjj 

- [0 N. K || dll +[ p C p ;; N.N.. dx dy] [ T j 1 = [0] 

(3.10) 

where the quantities enclosed within the square brackets 
indicate matrices . The closed boundary integral term of 
equation (3.10) can be substituted in terms of the 
boundary conditions of the problem noting that 

3T 

-K ax = q on dea ' t fl ux boundary (3.11a) 

= h(T-T f ) on convection (3.11b) 

boundary 

Now, substituting equations (3.11) into equation (3,10) 
the following equation is obtained 
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E 3N, 3N. 3N, 3N. . (e) 

K D f f I < T .. —H + jk — - -3 ! rfy H,, - ' r ™ 

e=1 l (e) lax Sy 9y - r y - - 3 


Eq 


Eh 


+ 2 0 N.q dl + E $ N. h ( N . T . ) dl 

e-1 (e) 1 e=1 1 0 ^ 


e= i 


Eh E 

” S 0 N h T 61 + 2 f p C {/ v N.N, dx dy^; r T ^ 

e=1 (e) e=1 ' p ^ e ' ) 1 3 o 


= r o 


(3.12) 


where , 


E = Number of internal or area elements 

Eq = Number of applied heat flux line elements 

Eh = Number of convective line elements . 

Rewriting the set of equation (3.12), 

2 lT^ e ^j + 1 q [Qi 6 ^] + 2 h f[H^ e ^ {T^ e H -[H^ 

;=1 e=1 " ° e=1 ^ 


+ z [m^] (t] -To] 

e=1 


(3.13) 


where , 


[K (e) ] 


, 3N. 3N. 3N. 3N, . 

// ( ... . — -H + — i. - y l dx dv 

^l3x 3x 3y 3y j y 


= Elemental conductivity matrix . 
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[Q B £) ]= (e) N i <1 « 

= Applied heat flux boundary vector* 

[H| e) ]= / h N N dl 

(e) x J 

= Left hand side matrix contribution from a 
convection boundary element 

C4 e) ] - , ' h T N dl 
(e) 

= Right hand side vector contribution from 
a convection boundary element. 

fM (e) ] = 1 ‘ N.N. dx dy 

(e) 1 J 

= Elemental mass matrix. 


In the present study, the elements used are eight 
noded isoparametric elements as shown in Fig. (3,1) whose 
boundaries are defined by parabolic curves in general . 

Since, the element sides can assume a curvilinear shape, 
complex boundaries of the solution domain can be accurately 
represented. Also, it is possible to define local coordinates 
X and V for an isoparametric element such that the global 
coordinates x,y of a point inside the element can be expressed 
in terms of the nodal coordinates and interpolation 
functions in a similar way as the temperature in equation 
(3.8). Therefore, 
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8 

x = E N. x. 


(3.14a) 


8 

y = 2 N y* (3.14b) 

i =1 A 


The interpolation functions I\L for the 8 -noded isoparametric 
element in terms of the local coordinates X and *), are 
given by the following relations, see ref. [9 1 . 

For corner nodes : 


N ± - l (1 


+ )(1+^ U) ( ^ + TJ i rj -1) 


(3.15 ) 


For mid -side nodes: 

\ = \ (1- ^ 2 )(1+ t ) ± V), if K. = 0 (3.16a) 


N ± = \ (1+ ^ HI-*? 2 ), if = 0 (3.16b) 

Using the interpolation functions of (3.16) in the finite 
element equations (3.13), one obtains the global matrices 
of the form : 

[M] fT} + {[KMRJ} [T] = {[H 2 ] -[Q B ]j (3.17) 

where , 

[ M ] is the Global mass matrix. 

[ k ] is the Global conductivity matrix. 



{ T J Is the vector containing time derivatives of 
the unknown nodal variable T. 

I 1 1 the vector containing unknown nodal 
variables of temperature. 

Equation (3.17) represents a set of ordinary differential 
equations of the variable T and its time derivatives . 
'-These sets of equations are solved by applying the finite 
difference technique. 


The steps involved are as follows 


For k time step equation (3.17) may be written as 
[«] (ij + {[ K] -[Hi! {T* ={[H 2 ] - [Q b 1 } (3.18) 


For the (k+1 b* 1 time step, it is 


. k+1 k+1 , k+1 

[M](T| + {[K] + [li,]} (T) ={[H 2 ] - [Q B ] j 

(3ll9) 

The finite difference scheme for the time derivative can 

now be written in a most general form as 

k+1 k+1 k . k 

(T } - ©• At { T } = \T ) + (1- 9 ) At {T J 


(3.20) 



■ 2 ~> 


where, 

is the required time step, 
k and k+1 refer to the k th and k+1 th time level. 

© is a constant, which decides the particular 
type of the scheme to be chosen, 

© = 0 gives the Explicit or Forward difference scheme. 

© =0.5 gives the Crank-Nicholson or Mid -difference 
scheme . 

© =1 .0 gives the Implicit or Backward difference 
scheme. 

In the present case, the implicit scheme is employed. The 
implicit method uses the time derivative evaluated at 
the (k+1) time level. Therefore, for each increment of 
time by At, the calculation of the nodal temperatures at 
the new time level requires the solution of a set of 
simultaneous algebraic equations. The solution is thus 
marched in time, in steps of At until the desired final time 
is reached. Although computations involving the implicit 
method are more complicated than that of the explicit method, 
the implicit formulation has the important advantage of 
being unconditionally stable [,10]. Larger values of At 
may therefore be used resulting in less computational time. 
Substituting equation (3.20) into equations (3.18) and (3.19) 
and rearranging, we get. 
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k+1 

( [M] + 9 At { [K ] +[ ft,] \ ) (t\ 

= ( [M]- (1- 9) At {[K] + [ft,] } ) (T} k 

+ At [H 2 ] - [Q b ]} (3.21) 

Since the vector (T} k is known for each k th time, the 
nodal temperature vector [T} k+1 can be calculated. 

3.5 CONVERSION OF THE PROBLEM TO AX I -SYMMETRIC FORMULATION 

In order to convert the problem to an axi -symmetric 
situation, the following modifications are carried out in 
the governing differential equation (3.5). 

Galerkin residue equation are : 

;; N ± ( v 2 T - ||) dv = [0] (3.22) 

where , 

dv = volume of a differential control volume under 
consideration. 

= 2nx dx dy for the axi -symmetric geometry. 

The equation (3.22) can be written as : 

JJ N ± { ( V 2 T - §|) } 2*x dx dy = [o 1 
or 

;/ ( «J. (N ± v T)- v N i . v T - N ± ) 2*x dx dy 

= [ o] 


(3.23) 



41 


The first term of equation (3.23) after the application of 
the divergence theorem, can be written as : 

// V . (N, y T) dv = 0 n. N. V T 2*x dl 

X — 

- // v N. . v T 2nx dx dy 

- l - 


Substituting the above expression into the equation (3.23) 
we get, 



3N . , 

yUdx dy] 



- [ / / N. N . x dx dy ] 




x dl 


0 j 


(3.24) 


The last term in equation (3.24) can be rewritten in terms 
of the applied boundary conditions of the problem. The 
resulting equations are of the form: 


E 

Z 

e=1 


3N. 

3N . 

■j-X + 

3x 

+ i q 

/ N 

e=1 

(e) ' 

E 


+ Z 

[ /> 

e=1 

(e) 


% 

3y 

q 


3} 

3y i 


x dl + 


N.N. x dx 

X J 


x dx dy] [ T 3 (e) ] 

^ I N.h (N.T ) x dl 
e=1 

dy][ Tj] -f 0] (3.25) 
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The set of simultaneous ordinary differential 
equations represented by (3.25), can be solved using the 
implicit finite difference procedure as discussed in the 
previous section. 

3.6 PHYSICAL REPRESENTATION OF THE PROBLEM : 

In order to proceed with the calculation for 
material removal, the geometry of the solution domain and 
the boundary conditions should be well defined. The 
physical domain, and the node numbering of the problem are 
done as shown in Fig. 3.2 and 3.3 respectively. The figures 
are self explanatory. 

3.7 DOMAIN VARIATION WITH TIME : 

Initially the domain is as shown in Fig. 3.3 and 
is considered to be at room temperature . After the 
discharge takes place for a given amount of time, the value 
of temperature at each node changes, depending on the type 
of boundary conditions, the domain is subjected to. Since 
we propose material removal to be primarily due to melting, 
the shape of the melting isotherm is drawn. A typical 
melting isotherm has a shape as shown in Fig. 3 A . The 
material above the isotherm (hatched area) is considered to 
be removed within the time step considered for computation. 
Now, for the next discharge, this melting isotherm becomes 



R Radius of spark 

Z Depth of penetration 


Fig. 3.2 Physical domain (work piece configuration) 
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1,2,3 - Node Numbers 

CD » (D » (J) ** Element Numbers 
T, 2, 3 - Boundary Node Numbers 


Fig. 3.3 Numbering of the grid of Fig. 3.2 
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Region above melting temperature or 
Material removed. 

Changing surfaces (boundaries) 

Melting isotherm 


Fig. 3.4 Melting isotherm and surface change 
representation. 
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"the surface of "the work piece on which “the sparking occurs , 
The nodes of the bottom surface (A3) are also moved a 
distance equal to those of the top surface, so that for any 
discharge the depth and radius of heat affected area for a 
spark remains relatively large (about 10 times) with respect 
to the spark radius . This process of updating the domain 
geometry is done everytime before applying the subsequent 
discharge. Thus, the FEM analysis is performed for the 
modified geometry after each time step. 



CHAPTER- IV 


RESULTS AND DISCUSSIONS : 

From the finite element formulation of the EDM 
process described in the previous chapter, results have 
been predicted for the temperature variation with respect 
to time and space in an H.S.S. work piece which is 
initially at room temperature. Assuming that the regions 
which attain a temperature above the upper limit of 
melting temperature (T^+A) are completely removed, the 
rate of material removal and the shape of the crater in 
the work piece for a typical spark have also been 
predicted . 

The results corresponding to a single spark have 
been obtained by considering total machining time of the 
order of typical pulse time with a very high power input. 

To get realistic results of an actual EDM problem a 
multispark model has also been solved considering 20 sparks 
with a continuous but low power input over a total 
machining time which is much larger than the pulse time. 

In Fig. 4.1 the amount of material removed by a 
single spark within a pulse time of 600 us has been 
compared with the analytical results of Snoeys and Dyck 
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as presented in Ref. [_ 5j . Various values of average 
discharge power during the pulse have been considered for 
comparison. Since the model presented [Ref . 5] takes into 
account partial evaporation and partial melting, the 
appropriate curve that needs to be compared with our 
result is the 100 */. melting curve. For the range of powers 
considered, the analytical and numerical results match 
approximately . The slight discrepancy observed between the 
results can be explained from the fact that no heat loss from 
the top surface is considered in the analytical result and 
also the heat flux is always applied at the top surface. 

The present formulation on the other hand accounts for 
heat loss from the top surface into the dielectric medium 
and also it updates the whole surface shape after each small 
time step At. 

The amount of material removed with time has been 
plotted for three different numerical cycle times in 
Fig. 4.2a. It is seen from the figure that the amount of 
material removed increases with time at a non-uniform rate . 

For certain times the. rate of machining is very high and for 
some others it is very low. This non uniformity in the rate of 
material removal can be attributed to the processes of 
temperature rise from room conditions upto the melting value 
and the phase change from solid to liquid, successively. 
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With respect to numerical cycle time it is observed that 
smaller the time step higher the value of material removed . 

When small time steps are used the molten material is 
removed more frequently and the newly formed surface is 
directly exposed to the heat source. This results in a higher 
overall material removal rate. In Fig, 4,2b the shape of 
the crater at the end of a total machining time of 2.5 sec. 
is shown for different numerical cycle times. For the reasons 
discussed above, the depth of penetration is more, for smaller 

At . 

In Fig. 4.3 the change in the shape of the crater with 
machining time for different time step At is presented. 

It is observed that in the initial stages of machining the 
crater extends deeper along the centre line of the spark. For 
later times, it spreads more in the radial direction upto 
spark radius.. For still larger machining time, the crater 
develops a taper along the sides. The width of the crater at 
the top surface of the workpiece is nearly equal to the spark 
diameter. In the present formulation a Gaussian distribution 
of heat flux intensity has been assumed within the spark. 

Due to this the rate of penetration along the centre line of 
the spark is higher for all times . It is also seen that for a 
smaller numerical cycle time the crater volume is larger at 
any instant of time . 
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In Fig. 4.4(a and. b) the effect of spark radius on 
material removal is shown. The plots indicate that at a 
fixed power level, the smaller the spark radius higher is the 
amount of material removed for a given time. For a smaller 
spark radius, the power intensity is higher and this leads 
to a higher rate of material removal , It is also seen that 
machining starts later if the spark size is larger. In 
Fig. 4.4b the shape of the crater after a machining time of 
2.5 sec. shows that for a smaller sized spark the width of 
the crater is less but the depth is much more than those of 
a larger sized spark. Thus it is clear that the radius of 
crater depends on the radius of spark but the depth of crater 
after a certain machining time depends on the average power 
intensity of the spark. These conclusions are further confirmed 
by the shapes of craters shown in Fig. 4.5 at various instances 
of time for two different spark radii. 

The effects of power per spark upon material removal 
is shown in Fig. 4.6. The chief influence of increasing the 
power per s£ark is that melting starts early . Infact there 
is an almost constant time gap between all the curves at all 
times. However, there is a very slight increase in slope with 
an increase in power per spark so that the material removal 
is slightly higher. The shapes of the craters at various 
instances of time, Fig. 4.7 , are shown for different power 
per spark. The crater shapes confirm the finding discussed 


above . 
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Radius i jjrn 



Fig. 4.5 Comparison of crater formation 
for different spark radii 
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In Fig. 4.8 "the variation of material removal rate 
(MRR) with power per spark has been presented, for a given 
spark radius and numerical cycle time. As expected MRR 
increases with power, since higher power input can supply 
more heat for melting the material. The MRR however, has a 
tendency to level off as power input per spark increases. This 
may be because of heat losses increasing at a fast rate . 

In Fig. 4.9 the effect of cycle time At upon MRR 
has been shown, with other parameters fixed. As cycle time 
decreases, MRR increases because of direct access of the work 
surface to the heat input more frequently . The MRR is seen to 
be very sensitive to changes in the value of the spark radius 
(Fig. 4.10). As the spark radius decreases the MRR increases 
very rapidly tending to approach the limit of a point source 
at small values of radius . Similar trends with variation in 
spark radius have been reported by Erden and Kaftanaglu [ "5 ]. 

The experimentally observed values of MRR by 
Indurkhya ( n) for three different machining conditions are 
also shown in the figure. It is seen that the range of 
variation of MRR is the same as in the present results .further 
validating our model . More detailed comparison of the experiment 
of Ref. [11] and present results is not possible due to the 
inavailability of important numerical parameters such as spark 
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Fig. 4.9 Variation of average MRR with ht 
(cycle time) 
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Fig. 4.10 Variation of average 


MRR with spark radius 
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radius and number of sparks in the gap for the experimental 
study . Also our theoretical study has a limitation of the 
continuous heat flux model which has been invoked to reduce 
the computational effort. The experimental machining conditions 
for which MRR values have been taken from Ref. [11] are 
summarized in Table 4.1 

The variation of material removal rate with machining 
time is shown in Fig. 4.11 for three different numerical cycle 
times . The MRR rapidly increases in the beginning after actual 
machining starts and a little while later it reaches a more or 
less steady value. There are always a few fluctuations in 
MRR for all time indicating the basic random nature of the 
EDM process . 

Fig. 4.12 gives the isotherms of different tempera- 
ture values. The plotted values show the movement of (melting) 
crater and it is fairly clear that the temperature increases 
very rapidly with depth in the (z) direction when compared 
to the radial (r) direction. 
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Temperature ( K ) 
a - 2000 

b - 1950 — Mel t i ng 

c - 1900 

d - 1875 

e - 1850 

f - 1800 


Fig. 4.12 Typical isotherms prevailing in 
the work piece 
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Table 1 Experimental MRR values for different 
machining conditions 


Curve 

ON time j 
JJS 

— „ . 

Tool radius 
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MRR | 

! 

mg /min j 
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1 

26.043 
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100 j 

6 

j 

24.230 1 

; c 

1 00 r 
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21.169 ! 
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5 . CONCLUSIONS AND SUGGESTIONS 

5.1 CONCLUSIONS : 

1 * The Physical formulation and the Finite Element 

formulation of the EDM process presented in this work 
appear to describe the actual phenomena reasonably 
accurately . The present analysis incorporates many 
improvements over previous theoretical works such as 
the accounting of the movement of the work surface, 
accounting of heat losses to the dielectric medium 
and variation of power intensity within the spark. 
Both single spark and multispark models have been 
developed to predict the rate of machining within a 
single pulse and also over many pulses upto a 
reasonable machining time. 

2. It is observed that the width of the spark depends 

on spark radius while the depth of penetration depends 
upon the power intensity. 

3. The numerical cycle time (time step, At) plays an 
important role in accurately predicting material 
removal . 

4. The material removal rate increases with power per 
spark, decreases with an increase in time step At. 
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It has a very strong dependence on spark radius 
appearing to reach the limit of infinity i.e. for a 
point source for very small spark radii. With respect 
to machining time the MRR increases rapidly immediately 
after the start of machining but soon it reaches a 
steady state value . 

5.2 SUGGESTIONS FOR FUTURE WORK : 

1 . The effect of evaporation of material can be 

considered by incorporating latent heat of boiling 
(evaporation) in addition to the latent heat of melting. 

2. The number of sparks can be assumed based on the 
surface roughness of tool and work piece materials. 

3. The heat source location can be assumed to be moving, 
with respect to time and space based on the relative 
surface roughness prevailing in the inter electrode 
gap. 

4. The temperature distribution in tool and dielectric 
medium can be predicted by extending the model to 
tool surface . Relative electrode wear can also be 
predicted. 
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